% ! problem10.m --> based on <<MATLAB codes for Finite Element Analysis>>
% ! chpater 7  using 2d-frame element

clear ;
clc
format default
p=struct();

% define model parameters

p.E = 210000;   % SI(mm)
p.A = 200;
p.P = [15000 10e6]';
p.Iz=2e8;

% define elems and nodes

p.nodes=[0 0 ;  % node's coordinates (x,y)
         0 6000 ;
         6000 6000 ;
         6000 0;];
p.elems=[1 2; 2 3; 3 4;];

p.node_num=size(p.nodes,1);
p.elem_num=size(p.elems,1);
p.node_Coord_x=p.nodes(:,1);
p.node_Coord_y=p.nodes(:,2);

% define problem's dimension = each node's dof  
p.problem_dimension=3;

% global degree of freedom number
p.global_dof_num=3*p.node_num;
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);

% define boundary conditions
% fixed dof 
p.fix_dof=[1 2 3  10 11 12]';
% load dof and amplitude
p.load_dof=[4 6]';

% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);

% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=p.problem_dimension*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 2d frame elem has 2 node,each node have 3 dof,
    % i-th node --> U_(3*i-2) and U_(3*i-1) and  U_(3*i) dof(displacement of node) 
    elem_dof=[connectivity(1)*3-2 connectivity(1)*3-1 connectivity(1)*3 connectivity(2)*3-2 connectivity(2)*3-1 connectivity(2)*3];

    node1_x=p.node_Coord_x(connectivity(1));
    node1_y=p.node_Coord_y(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2));
    node2_y=p.node_Coord_y(connectivity(2));
    
    elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2);
    l=(node2_x-node1_x)/elem_length;
    m=(node2_y-node1_y)/elem_length;

    % coordinate transform matrix: U_l=L*U_g
    L=[l m 0 0 0 0;
        -m l 0 0 0 0;
        0 0 1 0 0 0;
        0 0 0 l m 0;
        0 0 0 -m l 0;
        0 0 0 0 0 1;];

    % in local coordinates, the stiffness matrix of the frame element is obtained by com-
    % bination of the stiffness of the bar element and the Bernoulli beam element
    E=p.E;
    A=p.A;
    Iz=p.Iz;

    k_e_local=zeros(elem_global_dof_num,elem_global_dof_num);

    k_e_local(1,1)=(E*A)/elem_length;
    k_e_local(4,4)=(E*A)/elem_length;
    
    k_e_local(1,4)=-(E*A)/elem_length;
    k_e_local(4,1)=-(E*A)/elem_length;
    
    k_e_local(2,2)=12*E*Iz/(elem_length^3);
    k_e_local(2,3)=6*E*Iz/(elem_length^2);
    k_e_local(2,5)=12*E*Iz/(elem_length^3);
    k_e_local(2,6)=6*E*Iz/(elem_length^2);
    
    k_e_local(3,2)=6*E*Iz/(elem_length^2);
    k_e_local(5,2)=12*E*Iz/(elem_length^3);
    k_e_local(6,2)=6*E*Iz/(elem_length^2);
    
    k_e_local(3,3)=4*E*Iz/(elem_length);
    k_e_local(3,5)=-6*E*Iz/(elem_length*elem_length);
    k_e_local(3,6)=2*E*Iz/elem_length;
    
    k_e_local(6,3)=2*E*Iz/elem_length;
    k_e_local(5,3)=-6*E*Iz/(elem_length*elem_length);
    
    k_e_local(5,5)=12*E*Iz/(elem_length*elem_length*elem_length);
    k_e_local(5,6)=-6*E*Iz/(elem_length*elem_length);
    k_e_local(6,5)=-6*E*Iz/(elem_length*elem_length);
    
    k_e_local(6,6)=4*E*Iz/elem_length;

    k_e_global=L' *k_e_local *L;
    p.elemStiffs(:,:,i)=k_e_global;

    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e_global;

end

% apply boundary condition 
p.displacements(p.fix_dof)=0;
p.node_forces(p.load_dof)=p.P;

% solution KU=F
p = solutionStruct(p);

% show compute result
disp('displacements of dofs :')
disp([(1:p.global_dof_num)', p.displacements])

disp('forces of dofs :')
disp([(1:p.global_dof_num)',p.node_forces])

figure
axis equal

% plot undeformed mesh
draw_2Dtruss_mesh(p.nodes,p,'black')

deform_scale=500;
deformed_node_coordinate=[p.node_Coord_x+deform_scale*p.displacements(3*(1:p.node_num)-2), ...
                            p.node_Coord_y+deform_scale*p.displacements(3*(1:p.node_num)-1)];
% plot deformed mesh
draw_2Dtruss_mesh(deformed_node_coordinate,p,'red')

title('位移变形')
grid on

